function [zx] = interp3(z1, z2, z3, z4, x)
%
% Third Order Interplolation Function
%
%DESCRIPTION:
%This function computes the 3rd order interpolation to get the value zx at
%the position x with x€[0,1].
%
%PROTOTYPE:
% [zx] = interp3(z1, z2, z3, z4, x)
%
%--------------------------------------------------------------------------
% INPUTS:
%   z1        [1x1]       Node 1                     [-]
%   z2        [1x1]       Node 2                     [-]
%   z3        [1x1]       Node 3                     [-]
%   z4        [1x1]       Node 4                     [-]
%   x         [1x1]       Position                   [-]
%--------------------------------------------------------------------------
% OUTPUTS:
%   zx        [1x1]       Interpolated Value         [-]
%--------------------------------------------------------------------------
%
%NOTES:
% (none)
%
%CALLED FUNCTIONS:
% (none)
%
%UPDATES:
% (none)
%
%REFERENCES:
% [1] "Ionospheric Correction Algorithm for Galileo Single-Frequency Users"
%      - European GNSS (Galileo) Open Service
%
%AUTHOR(s):
%Luigi De Maria, Matteo D'Addazio, 2022
%

%% Main Code

if abs(2*x) < 1e-10
    zx = z2;
    
else
    delta = 2*x - 1;
    g1 = z3 + z2;
    g2 = z3 - z2;
    g3 = z4 + z1;
    g4 = (z4 - z1)/3;
    a0 = 9*g1 - g3;
    a1 = 9*g2 - g4;
    a2 = g3 - g1;
    a3 = g4 - g2;
    zx = 1/16 * (a0 + a1*delta + a2*delta^2 + a3*delta^3);
end

end